load packages

library(tidyverse)

source data

source("load_data.R")

standardize brain data

# roi distributions
dataset %>%
  filter(session == "all") %>%
  select(subjectID, con, condition, control, roi, process, meanPE) %>%
  unique() %>%
  filter(!is.na(roi)) %>%
  ggplot(aes(meanPE, fill = roi)) +
  geom_density(color = NA) +
  facet_wrap(~roi, ncol = 4) +
  theme_minimal() +
  theme(legend.position = "none")

# dot product distribution
dataset %>%
  filter(session == "all") %>%
  select(subjectID, con, condition, control, process, test, dotProduct) %>%
  unique() %>%
  ggplot(aes(dotProduct, fill = test)) +
  geom_density(color = NA) +
  facet_wrap(~process, ncol = 3, scales = "free") +
  theme_minimal()

# standardize
betas_std = betas %>%
  group_by(roi, session) %>%
  mutate(meanPE_std = scale(meanPE, center = TRUE, scale = TRUE),
         meanPE = ifelse(meanPE_std > 3 | meanPE_std < -3, NA, meanPE_std)) %>%
  ungroup()

dots_std = dots %>%
  group_by(map, test, mask, session) %>%
  mutate(dotProduct_std = scale(dotProduct, center = TRUE, scale = TRUE),
         dotProduct = ifelse(dotProduct_std > 3 | dotProduct_std < -3, 9999, dotProduct_std)) %>%
  ungroup()

# join data frames
dataset = full_join(betas_std, dots_std, by = c("subjectID", "con", "process", "condition", "control", "session")) %>%
  mutate(subjectID = as.character(subjectID)) %>%
  left_join(., ind_diffs, by = "subjectID") %>%
  ungroup() %>%
  mutate(subjectID = as.factor(subjectID),
         condition = as.factor(condition),
         control = as.factor(control),
         roi = as.factor(roi),
         process = as.factor(process),
         test = as.factor(test))

tidy data

model_data = dataset %>%
  filter(((!grepl("neuralsig", map) & test == "association" & mask == "masked") | (grepl("neuralsig", map) & mask == "unmasked")) & 
           session == "all" & control == "nature" & condition == "food") %>%
  unique() %>%
  group_by(process, subjectID) %>%
  mutate(meanProcessPEstd = mean(meanPE_std, na.rm = TRUE),
         processPE = paste0(process, "PE")) %>%
  select(-c(xyz, meanPE, meanPE_std, sdPE, dotProduct, map)) %>%
  spread(process, dotProduct_std) %>%
  spread(processPE, meanProcessPEstd) %>%
  group_by(subjectID) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  select(-c(roi, cravingPE, craving_regulationPE, test, age, gender, mask)) %>%
  unique()

run cross-validated models

BMI

data_ctrl = caret::trainControl(method = "cv", number = 5)

model_bmi_roi = caret::train(bmi ~ rewardPE + valuePE + cognitive_controlPE,
                     data = model_data,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

model_bmi_corr = caret::train(bmi ~ reward + value + cognitive_control,
                     data = model_data,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

model_bmi_cr = caret::train(bmi ~ craving_regulation,
                     data = model_data,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

model_bmi_combined = caret::train(bmi ~ rewardPE + valuePE + cognitive_controlPE + reward + value + cognitive_control + craving_regulation,
                     data = model_data,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

AIC(model_bmi_roi$finalModel, model_bmi_corr$finalModel, model_bmi_cr$finalModel, model_bmi_combined$finalModel)
summary(model_bmi_roi)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.2590 -2.1942 -0.2569  1.6289 14.9842 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          23.0416     0.2557  90.124   <2e-16 ***
## rewardPE             -1.0833     0.4651  -2.329   0.0207 *  
## valuePE               0.5376     0.3327   1.616   0.1073    
## cognitive_controlPE   0.1676     0.4604   0.364   0.7161    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.347 on 248 degrees of freedom
## Multiple R-squared:  0.02536,    Adjusted R-squared:  0.01357 
## F-statistic: 2.151 on 3 and 248 DF,  p-value: 0.09438
summary(model_bmi_corr)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.2667 -2.2289 -0.1992  1.5942 15.2710 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        22.9921     0.2712  84.790   <2e-16 ***
## reward             -0.7022     0.8992  -0.781    0.436    
## value               0.7088     0.7706   0.920    0.359    
## cognitive_control  -0.1118     0.4115  -0.272    0.786    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.378 on 248 degrees of freedom
## Multiple R-squared:  0.007218,   Adjusted R-squared:  -0.004791 
## F-statistic: 0.601 on 3 and 248 DF,  p-value: 0.6149
summary(model_bmi_cr)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.2408 -2.2708 -0.2915  1.6320 15.5800 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         22.9385     0.2555  89.782   <2e-16 ***
## craving_regulation   0.5522     0.3245   1.702     0.09 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.357 on 250 degrees of freedom
## Multiple R-squared:  0.01145,    Adjusted R-squared:  0.007499 
## F-statistic: 2.896 on 1 and 250 DF,  p-value: 0.09002
summary(model_bmi_combined)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.2589 -2.2038 -0.2826  1.8436 14.6980 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          22.7839     0.3155  72.226   <2e-16 ***
## rewardPE             -1.0780     0.5524  -1.951   0.0521 .  
## valuePE               0.2161     0.4468   0.484   0.6290    
## cognitive_controlPE   1.3243     1.0866   1.219   0.2241    
## reward               -0.4199     0.8981  -0.468   0.6405    
## value                 0.7507     0.8080   0.929   0.3538    
## cognitive_control    -0.8837     0.7803  -1.133   0.2585    
## craving_regulation    0.5476     0.3451   1.587   0.1139    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.339 on 244 degrees of freedom
## Multiple R-squared:  0.04526,    Adjusted R-squared:  0.01787 
## F-statistic: 1.652 on 7 and 244 DF,  p-value: 0.1215

body fat

data_ctrl = caret::trainControl(method = "cv", number = 5)

model_fat_roi = caret::train(fat ~ rewardPE + valuePE + cognitive_controlPE,
                     data = model_data,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

model_fat_corr = caret::train(fat ~ reward + value + cognitive_control,
                     data = model_data,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

model_fat_cr = caret::train(fat ~ craving_regulation,
                     data = model_data,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

model_fat_combined = caret::train(fat ~ rewardPE + valuePE + cognitive_controlPE + reward + value + cognitive_control + craving_regulation,
                     data = model_data,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

AIC(model_fat_roi$finalModel, model_fat_corr$finalModel, model_fat_cr$finalModel, model_fat_combined$finalModel)
summary(model_fat_roi)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.1941  -5.1750   0.9348   5.3876  19.5811 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          25.1642     0.6245  40.293   <2e-16 ***
## rewardPE             -1.5081     1.1370  -1.326    0.186    
## valuePE               0.6589     0.8121   0.811    0.418    
## cognitive_controlPE  -0.1986     1.1212  -0.177    0.860    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.104 on 245 degrees of freedom
## Multiple R-squared:  0.009479,   Adjusted R-squared:  -0.00265 
## F-statistic: 0.7816 on 3 and 245 DF,  p-value: 0.5052
summary(model_fat_corr)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.991  -5.246   1.268   5.267  18.942 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        24.9624     0.6573  37.976   <2e-16 ***
## reward             -2.5613     2.1623  -1.185    0.237    
## value               1.8902     1.8603   1.016    0.311    
## cognitive_control   0.1510     0.9927   0.152    0.879    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.101 on 245 degrees of freedom
## Multiple R-squared:  0.01016,    Adjusted R-squared:  -0.001957 
## F-statistic: 0.8385 on 3 and 245 DF,  p-value: 0.4739
summary(model_fat_cr)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.532  -5.337   1.246   4.997  19.444 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         25.3371     0.6214  40.771   <2e-16 ***
## craving_regulation   0.1318     0.7871   0.167    0.867    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.109 on 247 degrees of freedom
## Multiple R-squared:  0.0001135,  Adjusted R-squared:  -0.003935 
## F-statistic: 0.02804 on 1 and 247 DF,  p-value: 0.8671
summary(model_fat_combined)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.822  -5.097   1.056   5.663  18.591 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         24.82756    0.77580  32.003   <2e-16 ***
## rewardPE            -0.86398    1.35204  -0.639    0.523    
## valuePE              0.97665    1.09044   0.896    0.371    
## cognitive_controlPE  1.97281    2.65142   0.744    0.458    
## reward              -2.35581    2.19313  -1.074    0.284    
## value                1.52194    1.97758   0.770    0.442    
## cognitive_control   -0.87829    1.90988  -0.460    0.646    
## craving_regulation   0.07419    0.84436   0.088    0.930    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.134 on 241 degrees of freedom
## Multiple R-squared:  0.01823,    Adjusted R-squared:  -0.01028 
## F-statistic: 0.6394 on 7 and 241 DF,  p-value: 0.7231

SCA

tidy data

dots_sca = dots_std %>%
  filter((!grepl("neuralsig", map) &  mask == "masked") | (grepl("neuralsig", map) & mask == "unmasked")) %>%
  unite(variable, process, test, condition, control, sep = "_", remove = TRUE) %>%
  mutate(variable = sprintf("multivariate_%s", variable)) %>%
  filter(session == "all") %>%
  select(-c(map, con, mask, session, dotProduct)) %>%
  left_join(., ind_diffs) %>%
  select(-c(sample, DBIC_ID, age, gender))

dots_rest_assoc = dots_sca %>%
  filter(grepl("rest", variable) & grepl("snack|meal|dessert|food", variable) & grepl("association", variable)) %>%
  spread(variable, dotProduct_std) %>%
  unique() %>%
  ungroup()

dots_rest_uniform = dots_sca %>%
  filter(grepl("rest", variable) & grepl("snack|meal|dessert|food", variable) & grepl("uniformity", variable)) %>%
  spread(variable, dotProduct_std) %>%
  unique() %>%
  ungroup()

dots_nature_assoc = dots_sca %>%
  filter(grepl("nature", variable) & grepl("snack|meal|dessert|food", variable) & grepl("association", variable)) %>%
  spread(variable, dotProduct_std) %>%
  unique() %>%
  ungroup()

dots_nature_uniform = dots_sca %>%
  filter(grepl("nature", variable) & grepl("snack|meal|dessert|food", variable) & grepl("uniformity", variable)) %>%
  spread(variable, dotProduct_std) %>%
  unique() %>%
  ungroup()

betas_sca = betas_std %>%
  group_by(subjectID, process, condition, control) %>%
  mutate(meanProcessPEstd = mean(meanPE_std, na.rm = TRUE)) %>%
  unite(variable, process, condition, control, sep = "_", remove = TRUE) %>%
  mutate(variable = sprintf("univariate_%s", variable)) %>%
  filter(session == "all") %>%
  select(-c(con, session, xyz, roi, meanPE, sdPE, meanPE_std)) %>%
  left_join(., ind_diffs) %>%
  select(-c(sample, DBIC_ID, age, gender))

betas_rest = betas_sca %>%
  filter(grepl("rest", variable) & grepl("snack|meal|dessert|food", variable)) %>%
  unique() %>%
  spread(variable, meanProcessPEstd) %>%
  ungroup()

betas_nature = betas_sca %>%
  filter(grepl("nature", variable) & grepl("snack|meal|dessert|food", variable)) %>%
  unique() %>%
  spread(variable, meanProcessPEstd) %>%
  ungroup()

BMI

# set na.action for dredge
options(na.action = "na.fail")

betas

if (file.exists("bmi_sca_models_betas.RDS")) {
  betas_combined = readRDS("bmi_sca_models_betas.RDS")
} else {
  # betas > rest
  data = betas_rest %>%
    select(-fat) %>%
    na.omit()
  lm_predictors = paste(names(select(betas_rest, -c(subjectID, bmi, fat))), collapse = " + ")
  lm_formula = formula(paste0("bmi ~ ", lm_predictors, collapse = " + "))
  
  full_model = lm(lm_formula,
                  data = data)
  
  betas_rest_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
  
  # betas > nature
  data = betas_nature %>%
    select(-fat) %>%
    na.omit()
  lm_predictors = paste(names(select(betas_nature, -c(subjectID, bmi, fat))), collapse = " + ")
  lm_formula = formula(paste0("bmi ~ ", lm_predictors, collapse = " + "))
  
  full_model = lm(lm_formula,
                  data = data)
  
  betas_nature_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
  
  # combine models
  betas_combined = bind_rows(betas_rest_sca, betas_nature_sca) %>%
    select(AIC, delta, BIC, df, logLik, weight, `(Intercept)`, everything()) %>%
    arrange(AIC)
  
  # save
  saveRDS(betas_combined, "bmi_sca_models_betas.RDS")
}
# set AIC for null model you want to compare model AIC values to
null = betas_combined %>%
  arrange(df) %>%
  filter(df == 2)

# tidy for plotting
plot_data = betas_combined %>%
  filter(!(df == 2 & delta > .75)) %>%
  arrange(AIC) %>%
  #slice(1:10000) %>%
  mutate(specification = row_number(),
         better.fit = ifelse(AIC == null$AIC[1], "equal",
                      ifelse(AIC < null$AIC[1], "yes", "no"))) %>%
  gather(variable, val, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
  mutate(multivariate = ifelse(grepl("multivariate", variable) & !is.na(val), 1, NA),
         cognitive_control = ifelse(grepl("cognitive_control", variable) & !is.na(val), 1, NA),
         reward = ifelse(grepl("reward", variable) & !is.na(val), 1, NA),
         value = ifelse(grepl("value", variable) & !is.na(val), 1, NA),
         craving = ifelse(grepl("craving_a", variable) & !is.na(val), 1, NA),
         craving_regulation = ifelse(grepl("craving_regulation", variable) & !is.na(val), 1, NA),
         rest = ifelse(grepl("rest", variable) & !is.na(val), 1, NA),
         nature = ifelse(grepl("nature", variable) & !is.na(val), 1, NA),
         snack = ifelse(grepl("snack", variable) & !is.na(val), 1, NA),
         food = ifelse(grepl("food", variable) & !is.na(val), 1, NA),
         meal = ifelse(grepl("meal", variable) & !is.na(val), 1, NA),
         dessert = ifelse(grepl("dessert", variable) & !is.na(val), 1, NA),
         association_test = ifelse(grepl("association", variable) & !is.na(val), 1, NA),
         uniformity_test = ifelse(grepl("uniformity", variable) & !is.na(val), 1, NA)) %>%
  select(-c(variable, val)) %>%
  unique()

order = plot_data %>%
  arrange(AIC) %>%
  mutate(better.fit.num = ifelse(better.fit == "yes", 1, 0)) %>%
  gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, starts_with("better"))) %>%
  filter(!is.na(value)) %>%
  group_by(variable) %>%
  mutate(order = sum(better.fit.num)) %>%
  select(variable, order) %>%
  unique()

# variables included in model
variable_names = names(select(plot_data, -starts_with("better"), -specification, -AIC, -BIC, -df, -logLik, -delta, -weight))

# plot top panel
a = plot_data %>%
  ggplot(aes(specification, AIC, color = better.fit)) +
    geom_point(shape = "|", size = 4, alpha = .1) +
    geom_hline(yintercept = null$AIC, linetype = "dashed", color = "#5BBCD6") +
    scale_color_manual(values = c("black", "#F43C13")) + #"#5BBCD6", 
    scale_y_continuous(breaks = scales::pretty_breaks(4)) +
    labs(x = "", y = "AIC\n") +
    theme_minimal() +
    theme(legend.position = "none")

# plot bottom panel
plot_data_b = plot_data %>%
  gather(variable, val, eval(variable_names)) %>%
  left_join(., order, by = "variable") %>%
  mutate(variable = gsub("\\(Intercept\\)", "intercept", variable))

b = plot_data_b %>%
  ggplot(aes(specification, reorder(variable, order), color = better.fit)) +
    geom_point(data = subset(plot_data_b, !is.na(val)), alpha = .1, size = .25, position = position_jitter(width = 0.25, height = 0.25)) +
    scale_color_manual(values = c("black", "#F43C13")) + #"#5BBCD6", 
    labs(x = "\nspecification number", y = "variables\n") +
    theme_minimal() +
    theme(legend.position = "none")

cowplot::plot_grid(a, b, ncol = 1, align = "v")

dots

if (file.exists("bmi_sca_models_dots.RDS")) {
  dots_combined = readRDS("bmi_sca_models_dots.RDS")
} else {
  # dots association > rest
  data = dots_rest_assoc %>%
    select(-fat) %>%
    na.omit()
  lm_predictors = paste(names(select(dots_rest_assoc, -c(subjectID, bmi, fat))), collapse = " + ")
  lm_formula = formula(paste0("bmi ~ ", lm_predictors, collapse = " + "))
  
  full_model = lm(lm_formula,
                  data = data)
  
  clust = parallel::makeCluster(getOption("cl.cores", 7))
  invisible(parallel::clusterCall(clust, "library", character.only = TRUE))
  parallel::clusterExport(clust, "data")
  
  dots_rest_assoc_sca = MuMIn::pdredge(full_model, cluster = clust, rank = "AIC", extra = "BIC")
  parallel::stopCluster(clust)
  
  # dots uniformity > rest
  data = dots_rest_uniform %>%
    select(-fat) %>%
    na.omit()
  lm_predictors = paste(names(select(dots_rest_uniform, -c(subjectID, bmi, fat))), collapse = " + ")
  lm_formula = formula(paste0("bmi ~ ", lm_predictors, collapse = " + "))
  
  full_model = lm(lm_formula,
                  data = data)
  
  clust = parallel::makeCluster(getOption("cl.cores", 7))
  invisible(parallel::clusterCall(clust, "library", character.only = TRUE))
  parallel::clusterExport(clust, "data")
  
  dots_rest_uniform_sca = MuMIn::pdredge(full_model, cluster = clust, rank = "AIC", extra = "BIC")
  parallel::stopCluster(clust)
  
  # dots association > nature
  data = dots_nature_assoc %>%
    select(-fat) %>%
    na.omit()
  lm_predictors = paste(names(select(dots_nature_assoc, -c(subjectID, bmi, fat))), collapse = " + ")
  lm_formula = formula(paste0("bmi ~ ", lm_predictors, collapse = " + "))
  
  full_model = lm(lm_formula,
                  data = data)
  
  clust = parallel::makeCluster(getOption("cl.cores", 7))
  invisible(parallel::clusterCall(clust, "library", character.only = TRUE))
  parallel::clusterExport(clust, "data")
  
  dots_nature_assoc_sca = MuMIn::pdredge(full_model, cluster = clust, rank = "AIC", extra = "BIC")
  parallel::stopCluster(clust)
  
  # dots uniformity > nature
  data = dots_nature_uniform %>%
    select(-fat) %>%
    na.omit()
  lm_predictors = paste(names(select(dots_nature_uniform, -c(subjectID, bmi, fat))), collapse = " + ")
  lm_formula = formula(paste0("bmi ~ ", lm_predictors, collapse = " + "))
  
  full_model = lm(lm_formula,
                  data = data)
  
  clust = parallel::makeCluster(getOption("cl.cores", 7))
  invisible(parallel::clusterCall(clust, "library", character.only = TRUE))
  parallel::clusterExport(clust, "data")
  
  dots_nature_uniform_sca = MuMIn::pdredge(full_model, cluster = clust, rank = "AIC", extra = "BIC")
  parallel::stopCluster(clust)
  
  # combine models
  dots_combined = bind_rows(dots_rest_assoc_sca, dots_rest_uniform_sca, dots_nature_assoc_sca, dots_nature_uniform_sca) %>%
    select(AIC, delta, BIC, df, logLik, weight, `(Intercept)`, everything()) %>%
    arrange(AIC) 
  
  # save
  saveRDS(dots_combined, "bmi_sca_models_dots.RDS")
}

data_combined = bind_rows(betas_combined, dots_combined) %>%
  arrange(AIC)
data_combined = readRDS("bmi_sca_models.RDS")

combined plot

# set AIC for null model you want to compare model AIC values to
null = data_combined %>%
  arrange(df) %>%
  filter(df == 2)

# tidy for plotting
plot_data = data_combined %>%
  filter(!(df == 2 & delta > .75)) %>%
  arrange(AIC) %>%
  slice(1:10000) %>%
  mutate(specification = row_number(),
         better.fit = ifelse(AIC == null$AIC[1], "equal",
                      ifelse(AIC < null$AIC[1], "yes", "no"))) %>%
  gather(variable, val, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
  mutate(ROI = ifelse(grepl("univariate", variable) & !is.na(val), 1, NA),
         multivariate = ifelse(grepl("multivariate", variable) & !is.na(val), 1, NA),
         cognitive_control = ifelse(grepl("cognitive_control", variable) & !is.na(val), 1, NA),
         reward = ifelse(grepl("reward", variable) & !is.na(val), 1, NA),
         value = ifelse(grepl("value", variable) & !is.na(val), 1, NA),
         craving = ifelse(grepl("craving_a", variable) & !is.na(val), 1, NA),
         craving_regulation = ifelse(grepl("craving_regulation", variable) & !is.na(val), 1, NA),
         rest = ifelse(grepl("rest", variable) & !is.na(val), 1, NA),
         nature = ifelse(grepl("nature", variable) & !is.na(val), 1, NA),
         snack = ifelse(grepl("snack", variable) & !is.na(val), 1, NA),
         food = ifelse(grepl("food", variable) & !is.na(val), 1, NA),
         meal = ifelse(grepl("meal", variable) & !is.na(val), 1, NA),
         dessert = ifelse(grepl("dessert", variable) & !is.na(val), 1, NA),
         association_test = ifelse(grepl("association", variable) & !is.na(val), 1, NA),
         uniformity_test = ifelse(grepl("uniformity", variable) & !is.na(val), 1, NA)) %>%
  select(-c(variable, val)) %>%
  unique()

order = plot_data %>%
  arrange(AIC) %>%
  mutate(better.fit.num = ifelse(better.fit == "yes", 1, 0)) %>%
  gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, starts_with("better"))) %>%
  filter(!is.na(value)) %>%
  group_by(variable) %>%
  mutate(order = sum(better.fit.num)) %>%
  select(variable, order) %>%
  unique()

# variables included in model
variable_names = names(select(plot_data, -starts_with("better"), -specification, -AIC, -BIC, -df, -logLik, -delta, -weight))

# plot top panel
a = plot_data %>%
  ggplot(aes(specification, AIC, color = better.fit)) +
    geom_point(shape = "|", size = 4, alpha = .1) +
    geom_hline(yintercept = null$AIC, linetype = "dashed", color = "#5BBCD6") +
    scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
    scale_y_continuous(breaks = scales::pretty_breaks(4)) +
    labs(x = "", y = "AIC\n") +
    theme_minimal() +
    theme(legend.position = "none")

# plot bottom panel
b = plot_data %>%
  gather(variable, value, eval(variable_names)) %>%
  left_join(., order, by = "variable") %>%
  mutate(value = ifelse(!is.na(value), "|", ""),
         variable = gsub("\\(Intercept\\)", "intercept", variable)) %>%
  ggplot(aes(specification, reorder(variable, order), color = better.fit)) +
    geom_text(aes(label = value), alpha = .1) +
    scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
    labs(x = "\nspecification number", y = "variables\n") +
    theme_minimal() +
    theme(legend.position = "none")

p = cowplot::plot_grid(a, b, ncol = 1, align = "v")
ggsave(p, filename = "~/Desktop/plot.png", width = 25, height = 10)
# set AIC for null model you want to compare model AIC values to
null = data_combined %>%
  arrange(df) %>%
  filter(df == 2)

# tidy for plotting
plot_data = data_combined %>%
  filter(!(df == 2 & delta > .75)) %>%
  arrange(AIC) %>%
  slice(1:10000) %>%
  mutate(specification = row_number(),
         better.fit = ifelse(AIC == null$AIC[1], "equal",
                      ifelse(AIC < null$AIC[1], "yes", "no"))) %>%
  gather(variable, val, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
  mutate(ROI = ifelse(grepl("univariate", variable) & !is.na(val), 1, NA),
         multivariate = ifelse(grepl("multivariate", variable) & !is.na(val), 1, NA),
         cognitive_control = ifelse(grepl("cognitive_control", variable) & !is.na(val), 1, NA),
         reward = ifelse(grepl("reward", variable) & !is.na(val), 1, NA),
         value = ifelse(grepl("value", variable) & !is.na(val), 1, NA),
         craving = ifelse(grepl("craving_a", variable) & !is.na(val), 1, NA),
         craving_regulation = ifelse(grepl("craving_regulation", variable) & !is.na(val), 1, NA),
         rest = ifelse(grepl("rest", variable) & !is.na(val), 1, NA),
         nature = ifelse(grepl("nature", variable) & !is.na(val), 1, NA),
         snack = ifelse(grepl("snack", variable) & !is.na(val), 1, NA),
         food = ifelse(grepl("food", variable) & !is.na(val), 1, NA),
         meal = ifelse(grepl("meal", variable) & !is.na(val), 1, NA),
         dessert = ifelse(grepl("dessert", variable) & !is.na(val), 1, NA),
         association_test = ifelse(grepl("association", variable) & !is.na(val), 1, NA),
         uniformity_test = ifelse(grepl("uniformity", variable) & !is.na(val), 1, NA)) %>%
  select(-c(variable, val)) %>%
  unique()

order = plot_data %>%
  arrange(AIC) %>%
  mutate(better.fit.num = ifelse(better.fit == "yes", 1, 0)) %>%
  gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, starts_with("better"))) %>%
  filter(!is.na(value)) %>%
  group_by(variable) %>%
  mutate(order = sum(better.fit.num)) %>%
  select(variable, order) %>%
  unique()

# variables included in model
variable_names = names(select(plot_data, -starts_with("better"), -specification, -AIC, -BIC, -df, -logLik, -delta, -weight))

# plot top panel
a = plot_data %>%
  ggplot(aes(specification, AIC, color = better.fit)) +
    geom_point(shape = "|", size = 4, alpha = .1) +
    geom_hline(yintercept = null$AIC, linetype = "dashed", color = "#5BBCD6") +
    scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
    scale_y_continuous(breaks = scales::pretty_breaks(4)) +
    labs(x = "", y = "AIC\n") +
    theme_minimal() +
    theme(legend.position = "none")

# plot bottom panel
plot_data_b = plot_data %>%
  gather(variable, val, eval(variable_names)) %>%
  left_join(., order, by = "variable") %>%
  mutate(variable = gsub("\\(Intercept\\)", "intercept", variable))

b = plot_data_b %>%
  ggplot(aes(specification, reorder(variable, order), color = better.fit)) +
    geom_point(data = subset(plot_data_b, !is.na(val)), alpha = .1, size = .25, position = position_jitter(width = 0.25, height = 0.25)) +
    scale_color_manual(values = c("black", "#F43C13")) + #"#5BBCD6", 
    labs(x = "\nspecification number", y = "variables\n") +
    theme_minimal() +
    theme(legend.position = "none")

cowplot::plot_grid(a, b, ncol = 1, align = "v")

summarize number of better fitting models

# plot_data %>%
#   gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
#   mutate(better.fit = ifelse(better.fit == "yes", 1, 0),
#          var.better = ifelse(better.fit == 1 & !is.na(value), 1, 0),
#          variable = gsub("\\(Intercept\\)", "intercept", variable),
#          variable = gsub("health:", "health  x  ", variable),
#          variable = gsub("health", "food type", variable)) %>%
#   group_by(variable) %>%
#   summarize(sum.var = sum(var.better, na.rm = TRUE),
#             sum.all = sum(better.fit, na.rm = TRUE),
#             percent = (sum.var / sum.all) * 100) %>%
#   kable(format = "pandoc", digits = 2)

% fat

betas

if (file.exists("fat_sca_models_betas.RDS")) {
  betas_combined = readRDS("fat_sca_models_betas.RDS")
} else {
  # betas > rest
  data = betas_rest %>%
    select(-bmi) %>%
    na.omit()
  lm_predictors = paste(names(select(betas_rest, -c(subjectID, bmi, fat))), collapse = " + ")
  lm_formula = formula(paste0("fat ~ ", lm_predictors, collapse = " + "))
  
  full_model = lm(lm_formula,
                  data = data)
  
  betas_rest_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
  
  # betas > nature
  data = betas_nature %>%
    select(-bmi) %>%
    na.omit()
  lm_predictors = paste(names(select(betas_nature, -c(subjectID, bmi, fat))), collapse = " + ")
  lm_formula = formula(paste0("fat ~ ", lm_predictors, collapse = " + "))
  
  full_model = lm(lm_formula,
                  data = data)
  
  betas_nature_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
  
  # combine models
  betas_combined = bind_rows(betas_rest_sca, betas_nature_sca) %>%
    select(AIC, delta, BIC, df, logLik, weight, `(Intercept)`, everything()) %>%
    arrange(AIC)
  
  # save
  saveRDS(betas_combined, "fat_sca_models_betas.RDS")
}
# set AIC for null model you want to compare model AIC values to
null = betas_combined %>%
  arrange(df) %>%
  filter(df == 2)

# tidy for plotting
plot_data = betas_combined %>%
  filter(!(df == 2 & delta > .75)) %>%
  arrange(AIC) %>%
  #slice(1:10000) %>%
  mutate(specification = row_number(),
         better.fit = ifelse(AIC == null$AIC[1], "equal",
                      ifelse(AIC < null$AIC[1], "yes", "no"))) %>%
  gather(variable, val, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
  mutate(multivariate = ifelse(grepl("multivariate", variable) & !is.na(val), 1, NA),
         cognitive_control = ifelse(grepl("cognitive_control", variable) & !is.na(val), 1, NA),
         reward = ifelse(grepl("reward", variable) & !is.na(val), 1, NA),
         value = ifelse(grepl("value", variable) & !is.na(val), 1, NA),
         craving = ifelse(grepl("craving_a", variable) & !is.na(val), 1, NA),
         craving_regulation = ifelse(grepl("craving_regulation", variable) & !is.na(val), 1, NA),
         rest = ifelse(grepl("rest", variable) & !is.na(val), 1, NA),
         nature = ifelse(grepl("nature", variable) & !is.na(val), 1, NA),
         snack = ifelse(grepl("snack", variable) & !is.na(val), 1, NA),
         food = ifelse(grepl("food", variable) & !is.na(val), 1, NA),
         meal = ifelse(grepl("meal", variable) & !is.na(val), 1, NA),
         dessert = ifelse(grepl("dessert", variable) & !is.na(val), 1, NA),
         association_test = ifelse(grepl("association", variable) & !is.na(val), 1, NA),
         uniformity_test = ifelse(grepl("uniformity", variable) & !is.na(val), 1, NA)) %>%
  select(-c(variable, val)) %>%
  unique()

order = plot_data %>%
  arrange(AIC) %>%
  mutate(better.fit.num = ifelse(better.fit == "yes", 1, 0)) %>%
  gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, starts_with("better"))) %>%
  filter(!is.na(value)) %>%
  group_by(variable) %>%
  mutate(order = sum(better.fit.num)) %>%
  select(variable, order) %>%
  unique()

# variables included in model
variable_names = names(select(plot_data, -starts_with("better"), -specification, -AIC, -BIC, -df, -logLik, -delta, -weight))

# plot top panel
a = plot_data %>%
  ggplot(aes(specification, AIC, color = better.fit)) +
    geom_point(shape = "|", size = 4, alpha = .1) +
    geom_hline(yintercept = null$AIC, linetype = "dashed", color = "#5BBCD6") +
    scale_color_manual(values = c("black", "#F43C13")) + #"#5BBCD6", 
    scale_y_continuous(breaks = scales::pretty_breaks(4)) +
    labs(x = "", y = "AIC\n") +
    theme_minimal() +
    theme(legend.position = "none")

# plot bottom panel
plot_data_b = plot_data %>%
  gather(variable, val, eval(variable_names)) %>%
  left_join(., order, by = "variable") %>%
  mutate(variable = gsub("\\(Intercept\\)", "intercept", variable))

b = plot_data_b %>%
  ggplot(aes(specification, reorder(variable, order), color = better.fit)) +
    geom_point(data = subset(plot_data_b, !is.na(val)), alpha = .1, size = .25, position = position_jitter(width = 0.25, height = 0.25)) +
    scale_color_manual(values = c("black", "#F43C13")) + #"#5BBCD6", 
    labs(x = "\nspecification number", y = "variables\n") +
    theme_minimal() +
    theme(legend.position = "none")

cowplot::plot_grid(a, b, ncol = 1, align = "v")

dots

if (file.exists("fat_sca_models_dots.RDS")) {
  dots_combined = readRDS("fat_sca_models_dots.RDS")
} else {
  # dots association > rest
  data = dots_rest_assoc %>%
    select(-bmi) %>%
    na.omit()
  lm_predictors = paste(names(select(dots_rest_assoc, -c(subjectID, bmi, fat))), collapse = " + ")
  lm_formula = formula(paste0("fat ~ ", lm_predictors, collapse = " + "))
  
  full_model = lm(lm_formula,
                  data = data)
  
  dots_rest_assoc_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
  
  # dots uniformity > rest
  data = dots_rest_uniform %>%
    select(-bmi) %>%
    na.omit()
  lm_predictors = paste(names(select(dots_rest_uniform, -c(subjectID, bmi, fat))), collapse = " + ")
  lm_formula = formula(paste0("fat ~ ", lm_predictors, collapse = " + "))
  
  full_model = lm(lm_formula,
                  data = data)
  
  dots_rest_uniform_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
  
  # dots association > nature
  data = dots_nature_assoc %>%
    select(-bmi) %>%
    na.omit()
  lm_predictors = paste(names(select(dots_nature_assoc, -c(subjectID, bmi, fat))), collapse = " + ")
  lm_formula = formula(paste0("fat ~ ", lm_predictors, collapse = " + "))
  
  full_model = lm(lm_formula,
                  data = data)
  
  dots_nature_assoc_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
  
  # dots association > nature
  data = dots_nature_uniform %>%
    select(-bmi) %>%
    na.omit()
  lm_predictors = paste(names(select(dots_nature_uniform, -c(subjectID, bmi, fat))), collapse = " + ")
  lm_formula = formula(paste0("fat ~ ", lm_predictors, collapse = " + "))
  
  full_model = lm(lm_formula,
                  data = data)
  
  dots_nature_uniform_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
  
  # combine models
  dots_combined = bind_rows(dots_rest_assoc_sca, dots_rest_uniform_sca, dots_nature_assoc_sca, dots_nature_uniform_sca) %>%
    select(AIC, delta, BIC, df, logLik, weight, `(Intercept)`, everything()) %>%
    arrange(AIC) 
  
  # save
  saveRDS(dots_combined, "fat_sca_models_dots.RDS")
}

data_combined = bind_rows(betas_combined, dots_combined) %>%
    arrange(AIC)
# set AIC for null model you want to compare model AIC values to
null = dots_combined %>%
  arrange(df) %>%
  filter(df == 2)

# tidy for plotting
plot_data = dots_combined %>%
  filter(!(df == 2 & delta > .75)) %>%
  arrange(AIC) %>%
  slice(1:10000) %>%
  mutate(specification = row_number(),
         better.fit = ifelse(AIC == null$AIC[1], "equal",
                      ifelse(AIC < null$AIC[1], "yes", "no"))) %>%
  gather(variable, val, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
  mutate(ROI = ifelse(grepl("univariate", variable) & !is.na(val), 1, NA),
         multivariate = ifelse(grepl("multivariate", variable) & !is.na(val), 1, NA),
         cognitive_control = ifelse(grepl("cognitive_control", variable) & !is.na(val), 1, NA),
         reward = ifelse(grepl("reward", variable) & !is.na(val), 1, NA),
         value = ifelse(grepl("value", variable) & !is.na(val), 1, NA),
         craving = ifelse(grepl("craving_a", variable) & !is.na(val), 1, NA),
         craving_regulation = ifelse(grepl("craving_regulation", variable) & !is.na(val), 1, NA),
         rest = ifelse(grepl("rest", variable) & !is.na(val), 1, NA),
         nature = ifelse(grepl("nature", variable) & !is.na(val), 1, NA),
         snack = ifelse(grepl("snack", variable) & !is.na(val), 1, NA),
         food = ifelse(grepl("food", variable) & !is.na(val), 1, NA),
         meal = ifelse(grepl("meal", variable) & !is.na(val), 1, NA),
         dessert = ifelse(grepl("dessert", variable) & !is.na(val), 1, NA),
         association_test = ifelse(grepl("association", variable) & !is.na(val), 1, NA),
         uniformity_test = ifelse(grepl("uniformity", variable) & !is.na(val), 1, NA)) %>%
  select(-c(variable, val)) %>%
  unique()

order = plot_data %>%
  arrange(AIC) %>%
  mutate(better.fit.num = ifelse(better.fit == "yes", 1, 0)) %>%
  gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, starts_with("better"))) %>%
  filter(!is.na(value)) %>%
  group_by(variable) %>%
  mutate(order = sum(better.fit.num)) %>%
  select(variable, order) %>%
  unique()

# variables included in model
variable_names = names(select(plot_data, -starts_with("better"), -specification, -AIC, -BIC, -df, -logLik, -delta, -weight))

# plot top panel
a = plot_data %>%
  ggplot(aes(specification, AIC, color = better.fit)) +
    geom_point(shape = "|", size = 4, alpha = .1) +
    geom_hline(yintercept = null$AIC, linetype = "dashed", color = "#5BBCD6") +
    scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
    scale_y_continuous(breaks = scales::pretty_breaks(4)) +
    labs(x = "", y = "AIC\n") +
    theme_minimal() +
    theme(legend.position = "none")

# plot bottom panel
plot_data_b = plot_data %>%
  gather(variable, val, eval(variable_names)) %>%
  left_join(., order, by = "variable") %>%
  mutate(variable = gsub("\\(Intercept\\)", "intercept", variable))

b = plot_data_b %>%
  ggplot(aes(specification, reorder(variable, order), color = better.fit)) +
    geom_point(data = subset(plot_data_b, !is.na(val)), alpha = .1, size = .25, position = position_jitter(width = 0.25, height = 0.25)) +
    scale_color_manual(values = c("black", "#F43C13")) + #"#5BBCD6", 
    labs(x = "\nspecification number", y = "variables\n") +
    theme_minimal() +
    theme(legend.position = "none")

cowplot::plot_grid(a, b, ncol = 1, align = "v")

combined plot

# set AIC for null model you want to compare model AIC values to
null = data_combined %>%
  arrange(df) %>%
  filter(df == 2)

# tidy for plotting
plot_data = data_combined %>%
  filter(!(df == 2 & delta > .75)) %>%
  arrange(AIC) %>%
  slice(1:10000) %>%
  mutate(specification = row_number(),
         better.fit = ifelse(AIC == null$AIC[1], "equal",
                      ifelse(AIC < null$AIC[1], "yes", "no"))) %>%
  gather(variable, val, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
  mutate(ROI = ifelse(grepl("univariate", variable) & !is.na(val), 1, NA),
         multivariate = ifelse(grepl("multivariate", variable) & !is.na(val), 1, NA),
         cognitive_control = ifelse(grepl("cognitive_control", variable) & !is.na(val), 1, NA),
         reward = ifelse(grepl("reward", variable) & !is.na(val), 1, NA),
         value = ifelse(grepl("value", variable) & !is.na(val), 1, NA),
         craving = ifelse(grepl("craving_a", variable) & !is.na(val), 1, NA),
         craving_regulation = ifelse(grepl("craving_regulation", variable) & !is.na(val), 1, NA),
         rest = ifelse(grepl("rest", variable) & !is.na(val), 1, NA),
         nature = ifelse(grepl("nature", variable) & !is.na(val), 1, NA),
         snack = ifelse(grepl("snack", variable) & !is.na(val), 1, NA),
         food = ifelse(grepl("food", variable) & !is.na(val), 1, NA),
         meal = ifelse(grepl("meal", variable) & !is.na(val), 1, NA),
         dessert = ifelse(grepl("dessert", variable) & !is.na(val), 1, NA),
         association_test = ifelse(grepl("association", variable) & !is.na(val), 1, NA),
         uniformity_test = ifelse(grepl("uniformity", variable) & !is.na(val), 1, NA)) %>%
  select(-c(variable, val)) %>%
  unique()

order = plot_data %>%
  arrange(AIC) %>%
  mutate(better.fit.num = ifelse(better.fit == "yes", 1, 0)) %>%
  gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, starts_with("better"))) %>%
  filter(!is.na(value)) %>%
  group_by(variable) %>%
  mutate(order = sum(better.fit.num)) %>%
  select(variable, order) %>%
  unique()

# variables included in model
variable_names = names(select(plot_data, -starts_with("better"), -specification, -AIC, -BIC, -df, -logLik, -delta, -weight))

# plot top panel
a = plot_data %>%
  ggplot(aes(specification, AIC, color = better.fit)) +
    geom_point(shape = "|", size = 4, alpha = .1) +
    geom_hline(yintercept = null$AIC, linetype = "dashed", color = "#5BBCD6") +
    scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
    scale_y_continuous(breaks = scales::pretty_breaks(4)) +
    labs(x = "", y = "AIC\n") +
    theme_minimal() +
    theme(legend.position = "none")

# plot bottom panel
plot_data_b = plot_data %>%
  gather(variable, val, eval(variable_names)) %>%
  left_join(., order, by = "variable") %>%
  mutate(variable = gsub("\\(Intercept\\)", "intercept", variable))

b = plot_data_b %>%
  ggplot(aes(specification, reorder(variable, order), color = better.fit)) +
    geom_point(data = subset(plot_data_b, !is.na(val)), alpha = .1, size = .25, position = position_jitter(width = 0.25, height = 0.25)) +
    scale_color_manual(values = c("black", "#F43C13")) + #"#5BBCD6", 
    labs(x = "\nspecification number", y = "variables\n") +
    theme_minimal() +
    theme(legend.position = "none")

cowplot::plot_grid(a, b, ncol = 1, align = "v")

summarize number of better fitting models

# plot_data %>%
#   gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
#   mutate(better.fit = ifelse(better.fit == "yes", 1, 0),
#          var.better = ifelse(better.fit == 1 & !is.na(value), 1, 0),
#          variable = gsub("\\(Intercept\\)", "intercept", variable),
#          variable = gsub("health:", "health  x  ", variable),
#          variable = gsub("health", "food type", variable)) %>%
#   group_by(variable) %>%
#   summarize(sum.var = sum(var.better, na.rm = TRUE),
#             sum.all = sum(better.fit, na.rm = TRUE),
#             percent = (sum.var / sum.all) * 100) %>%
#   kable(format = "pandoc", digits = 2)